function f = npqiv_ann_log(guess,Y,X1,X2,P,kh,tau,approx,k)

% Compute the criterion function for each parameter to estimate
% the SPQIV model using Chen and Pouzo's (2011) PSMD estimator.
% Use Gaussian Radial Basis ANN(kh) to approximate the function of interest.
% Written by Eduardo Souza-Rodrigues.
%
% INPUTS:   a) guess   = [alpha; mu; sigma; beta]
%           b) alpha   = vector of weights for h(w)
%           c) mu      = vector of means fo gaussian terms
%           d) sigma   = vector of stdvs fo gaussian terms
%           c) beta    = vector of finite dimensional parameters
%           d) Y       = endogenous variable (T x 1) vector
%           e) X1      = special regressor (T x 1) matrix
%           f) X2      = other endogenous regressors (T x dx) matrix 
%                        (do not include a constant term)
%           g) P       = matrix of instruments (T x J) matrix
%           h) kh      = dimension of Gaussian Radial Basis ANN(kh) + 1
%           i) lambda1 = penalization parameter
%           j) lambda2 = penalization parameter
%           k) tau     = estimate at tau-quantile
%           l) approx  = 0 use indicator function for residual function
%                      = 1 use smooth version for residual function
%           m) k       = smooth parameter for 'approx' = 1.
%
% OUTPUT: f is the value of the Criterion Function


% Important values
T    = size(Y,1);
dim  = size(guess,1);
   
% Separate [alpha, mu, sigma, beta] in guess
alpha = guess(1:kh,1);
mu    = guess(kh+1:2*kh,1);
sigma = guess(2*kh+1:3*kh,1);
beta  = guess(3*kh+1:dim,1);

% Compute single-index W
W = X2 * beta - X1;

% Basis function
G    = zeros(T,kh);
for j=1:kh
    G(:,j) = normcdf((Y-mu(j,1)*ones(T,1))/sigma(j,1));
end

% Use indicator function or a smoothed version for the residual function
if approx == 0      % Indicator
    R = (G*alpha-W<=0)-tau*ones(T,1);
elseif approx == 1  % Smooth version
    R = 0.5 + 0.5*(-tanh(k*(G*alpha-W)))-tau*ones(T,1);
end
% Approximation for indicator (x<0) is H(x) = 0.5 + 0.5*(-tanh(k*x)),
% where k is the smoothing factor (the smaller, the smoothier)

% Compute approximation for m(Z) using P
m = P*((P'*P)\(P'*R));

% Compute average of m squared (using identity matrix as weighting matrix)
M = mean(m.^2);

% Compute the Penalization Terms:
% (a) Smoothness
% Evaluation points
%evw = nodeunif(T,lb,ub);
% D = matrix with second derivative of B(.) at the evaluating points
%D = splibas(knots,0,kh,evw,2);
% C = integral over min(W) and max(W) of the terms D*D'
% Using the empirical measure for the integral, C becomes equivalent to (D'*D)/T
%C = (D'*D)/T;

% (b) Violation of monotonicity: h'>0
%Dh    = zeros(T,1);
%for j=1:kh
%    Dtemp = (alpha(j+1,1)/sigma(j,1))*normpdf((W-mu(j,1)*ones(T,1))/sigma(j,1));
%    Dh = Dh + Dtemp;
%end
%V3 = max(-Dh,0);
%E3 = mean(V3);

% Compute the Criterion Function:
f = M;
%f = M + lambda1 * alpha'*C*alpha + lambda2 * E;
%f = M + lambda2 * E2;
%f = M + lambda2 * E2 + lambda3 * E3;
